function bc = savecon(b, w)
% savecon  - save contrast computation taking missing values into account
%
% FORMAT:       bc = savecon(b, w)
%
% Input fields:
%
%       b           beta estimates
%       w           contrast weighting vector
%
% Output fields:
%
%       bc          beta contrast

% Version:  v0.7g
% Build:    9040222
% Date:     Apr-02 2009, 10:18 PM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
   ~isnumeric(b) || ...
   ~isa(w, 'double') || ...
    size(b, ndims(b)) ~= numel(w) || ...
    any(isinf(w(:)) || isnan(w(:)))
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument.' ...
    );
end

% find weights that are ~= 0
w = w(:)';
wn = find(w < 0);
wp = find(w > 0);

% check contrast
if ~isempty(wn) && ...
   ~isempty(wp)
    if sum(w([wn, wp])) ~= 0
        warning( ...
            'BVQXtools:Warning', ...
            'Differential contrast sum ~= 0!' ...
        );
    end
end

% re-weigh both parts to one
w(wn) = (-1 / sum(w(wn))) .* w(wn);
w(wp) = (1 / sum(w(wp))) .* w(wp);

% build contrast output
sb = size(b);
sb(end) = 1;
bn = zeros(sb);
bp = zeros(sb);
bs = zeros(sb);

% build subsref argument
sr = struct('type', '()', 'subs', {repmat({':'}, 1, numel(sb))});

% start with negative weights
for wc = wn
    sr.subs{end} = wc;
    bt = double(subsref(b, sr));
    bt(isinf(bt) | isnan(bt)) = 0;
    bn = bn + w(wc) .* bt;
    bs = bs + w(wc) .* (bt ~= 0);
end
bn = bn ./ abs(bs);
bn(bs == 0) = 0;

% continue with positive weights
bs = zeros(sb);
for wc = wp
    sr.subs{end} = wc;
    bt = double(subsref(b, sr));
    bt(isinf(bt) | isnan(bt)) = 0;
    bp = bp + w(wc) .* bt;
    bs = bs + w(wc) .* (bt ~= 0);
end
bp = bp ./ abs(bs);
bp(bs == 0) = 0;

% combine
bc = bn + bp;
